
#############################################################
####  Replication Code for Carney and Kaslovsky (2016)   ####
####                 April, 26, 2016                     ####
#############################################################


library(foreign)
library(readstata13)
library(lme4)
library(plyr)
library(stargazer)
library(mvtnorm)

setwd("~/downloads")
mydata <- read.dta("jop data-table 1 (2).dta")
ADA2008 <- read.dta13("Replication Data with ADA .dta")

###################################################
#        Recreating table 1 of Ellis (2013)       #
###################################################

## Step 1. Create unique district ID
mydata$unique.district <- paste(mydata$State_Code__v206_, mydata$CD_Code__v250_, sep=" - ")

## Step 2. Run a Mixed Model 
mle <- lmer(distancebasic ~ incomecenter + black + latino + votereported + actcenter + knowcenter + inparty + (1| unique.district), data = mydata)
summary(mle)

# Need to change the interaction variable names so that the code runs
mydata2 <- rename(mydata, c("_IinXun"="IinXun", "_IinXab"="IinXab", "_IinXgi"="IinXgi", "_IinXhh"="IinXhh", "_IinXre"="IinXre"))


mle2 <- lmer(distancebasic ~ incomecenter + black + latino + votereported + actcenter + knowcenter + inparty + I(hhincenter/10000) + ginicenter + rephouse + unioncenter + abspvicenter + I(incomecenter*(hhincenter/10000)) + IinXgi +  IinXab + IinXun  + IinXre + (1| unique.district), data = mydata2)
summary(mle2)

names(mydata)

mle3 <- lmer(keyvote ~ incomecenter + black + latino + votereported + actcenter + knowcenter + inparty + (1| unique.district), data = mydata)
summary(mle3)

mle4 <- lmer(keyvote ~ incomecenter + black + latino + votereported + actcenter + knowcenter + inparty + I(hhincenter/10000) + ginicenter + rephouse + unioncenter + abspvicenter + I(incomecenter*(hhincenter/10000))  + IinXgi +  IinXab + IinXun  + IinXre + (1| unique.district), data = mydata2)
summary(mle4)

# Put it in a table
x <- stargazer(mle,mle2,mle3,mle4, title="Replication of Table 1")


#################################################
#    Rerunning analysis with ADA Scores        #
################################################

#use new ADA dataset: ADA2008

##Rescaling the ADA Scores so that liberal is a lower score and conservative is higher from 0 to 100
ADA2008$RescaledADA <- rep(NA,nrow(ADA2008))
ADA2008$RescaledADA <- 100 - ADA2008$rating 

ADA2008$unique.district <- paste(ADA2008$State_Code__v206_, ADA2008$CD_Code__v250_, sep=" - ")
ADA2008$distancebasic2 <- abs(ADA2008$RescaledADA - ADA2008$ideol100pt)

library(stargazer)
library(lme4)
mle <- lmer(distancebasic2 ~ incomecenter + black + latino + votereported + actcenter + knowcenter + inparty + (1| unique.district), data = ADA2008)
summary(mle)


# Need to change the interaction variable names so that the code runs
library(plyr)
mydata2 <- rename(ADA2008, c("_IinXun"="IinXun", "_IinXab"="IinXab", "_IinXgi"="IinXgi", "_IinXhh"="IinXhh", "_IinXre"="IinXre"))


mle2 <- lmer(distancebasic2 ~ incomecenter + black + latino + votereported + actcenter + knowcenter + inparty + I(hhincenter/10000) + ginicenter + rephouse + unioncenter + abspvicenter + I(incomecenter*(hhincenter/10000)) + IinXgi +  IinXab + IinXun  + IinXre + (1| unique.district), data = mydata2)
summary(mle2)

# Put it in a table
library(stargazer)
x <- stargazer(mle,mle2, title="Table 2")


##########                                             #################
#       Ideological Distance: median household income                  #
##########                                            ##################

par(mfrow=c(1,2))  #put the graphs next to each other 


beta1 <- -0.244615 
beta3 <- -0.066706 # setting the coefficients 

margEff <- function(z){
  return(beta1 + z*beta3)
}

beta3
beta1

mine<- na.omit(mydata2$hhincenter/10000)


### Figure out the lower and upper x-axis bounds for plotting
xlb <- min(mine)
xup <- max(mine)

### Create a blank plot to add your lines to
plot(x=c(), y=c(), ylim=c(-1, 1), xlim=c(xlb, xup), 
     xlab="Median Household Income", ylab="Marginal Effect of Constituent Income", 
     main="")

### Add a horizontal line for ``zero"
abline(h=0, lty=3)

### Add your lines using the ``curve" function
curve(margEff, from = xlb, to=xup, lwd=2, add=T)

##### Calculate confidence intervals
#### Get the variance-covariance matrix
covmat <- vcov(mle2)
covmat

#### Get variances/covariances
vcovBeta1 <- covmat["incomecenter","incomecenter"]
vcovBeta3 <- covmat["I(incomecenter * (hhincenter/10000))","I(incomecenter * (hhincenter/10000))"]
covBeta1Beta3 <- covmat["incomecenter","I(incomecenter * (hhincenter/10000))"]

### Make a function that calculates variance of marginal effect given Z
margEffVar <- function(z){
  return(vcovBeta1 + (z^2)*vcovBeta3 + 2*z*covBeta1Beta3)
}
### Make a function that calculates the lower bound of the 95% confidence interval
margEff_95CI_LB <- function(z){
  return(margEff(z) - qnorm(.975)*sqrt(margEffVar(z)))
}
### Make a function that calculates the upper bound of the 95% confidence interval
margEff_95CI_UB <- function(z){
  return(margEff(z) + qnorm(.975)*sqrt(margEffVar(z)))
}

##### Add the upper and lower CI curves to your plot
# Lower
curve(margEff_95CI_LB, from=xlb, to=xup, lty=2, add=T)
# Upper
curve(margEff_95CI_UB, from=xlb, to=xup, lty=2, add=T)

#### Add a vertical line for the mean of the moderating variable
abline(v=mean(mine), lwd=2, lty=2, col="red")

#### Add a rug plot
rug(mine)

#### Make a transparent histogram
#### This is a function that turns any color semi-transparent
#### ``alpha" controls the degree of transparency
makeTransparent<-function(someColor, alpha=100){
  newColor<-col2rgb(someColor)
  apply(newColor, 2, function(curcoldata){rgb(red=curcoldata[1], green=curcoldata[2],
                                              blue=curcoldata[3],alpha=alpha, maxColorValue=255)})
}

#### Make a transparent grey color
hist_col = makeTransparent("grey")

# Tell R to add any new plots to the current one
par(new=T)
# Add a histogram with no axes and with the transparent color ``hist_col"
hist(mine, axes=F, xlab="", ylab="",main="", border=hist_col, col=hist_col)
# Tell R to stop adding new plots to the current one
par(new=F)


#################################################################
#       Ideological Distance:  District Competitiveness         #
#################################################################

mle2 <- lmer(distancebasic ~ incomecenter + black + latino + votereported + actcenter + knowcenter + inparty + I(hhincenter/10000) + ginicenter + rephouse + unioncenter + abspvicenter + I(incomecenter*(hhincenter/10000)) + ginicenter*incomecenter +  abspvicenter*incomecenter + unioncenter*incomecenter  + rephouse*incomecenter + (1| unique.district), data = mydata2)
summary(mle2)

beta1 <- -0.244615 
beta3 <- -0.015288 

margEff <- function(z){
  return(beta1 + z*beta3)
}

mine<- na.omit(mydata2$abspvicenter)


### Figure out the lower and upper x-axis bounds for plotting
xlb <- min(mine)
xup <- max(mine)

### Create a blank plot to add your lines to
plot(x=c(), y=c(), ylim=c(-1, 1), xlim=c(xlb, xup), 
     xlab="District Competitiveness", ylab="Marginal Effect of Constituent Income", 
     main="")

### Add a horizontal line for ``zero"
abline(h=0, lty=3)

### Add your lines using the ``curve" function
curve(margEff, from = xlb, to=xup, lwd=2, add=T)

##### Calculate confidence intervals
#### Get the variance-covariance matrix
covmat <- vcov(mle2)
covmat

#### Get variances/covariances
vcovBeta1 <- covmat["incomecenter","incomecenter"]
vcovBeta3 <- covmat["incomecenter:abspvicenter","incomecenter:abspvicenter"]
covBeta1Beta3 <- covmat["incomecenter","incomecenter:abspvicenter"]

### Make a function that calculates variance of marginal effect given Z
margEffVar <- function(z){
  return(vcovBeta1 + (z^2)*vcovBeta3 + 2*z*covBeta1Beta3)
}
### Make a function that calculates the lower bound of the 95% confidence interval
margEff_95CI_LB <- function(z){
  return(margEff(z) - qnorm(.975)*sqrt(margEffVar(z)))
}
### Make a function that calculates the upper bound of the 95% confidence interval
margEff_95CI_UB <- function(z){
  return(margEff(z) + qnorm(.975)*sqrt(margEffVar(z)))
}

##### Add the upper and lower CI curves to your plot
# Lower
curve(margEff_95CI_LB, from=xlb, to=xup, lty=2, add=T)
# Upper
curve(margEff_95CI_UB, from=xlb, to=xup, lty=2, add=T)

#### Add a vertical line for the mean of the moderating variable
abline(v=mean(mine), lwd=2, lty=2, col="red")

#### Add a rug plot
rug(mine)

#### Make a transparent histogram
#### This is a function that turns any color semi-transparent
#### ``alpha" controls the degree of transparency
makeTransparent<-function(someColor, alpha=100){
  newColor<-col2rgb(someColor)
  apply(newColor, 2, function(curcoldata){rgb(red=curcoldata[1], green=curcoldata[2],
                                              blue=curcoldata[3],alpha=alpha, maxColorValue=255)})
}

#### Make a transparent grey color
hist_col = makeTransparent("grey")

par(new=T)

# Add a histogram with no axes and with the transparent color ``hist_col"
hist(mine, axes=F, xlab="", ylab="",main="", border=hist_col, col=hist_col)
# Tell R to stop adding new plots to the current one
par(new=F)



############################################
#        Key Vote Variable Plots           #
############################################


#####################################################
#        Key Vote: median household income          #
#####################################################

par(mfrow=c(1,2))

beta1 <- -0.344473
beta3 <- 0.077905 

margEff <- function(z){
  return(beta1 + z*beta3)
}

beta3
beta1

mine<- na.omit(mydata2$hhincenter/10000)


### Figure out the lower and upper x-axis bounds for plotting
xlb <- min(mine)
xup <- max(mine)

### Create a blank plot to add your lines to
plot(x=c(), y=c(), ylim=c(-1, 1), xlim=c(xlb, xup), 
     xlab="Median Household Income", ylab="Marginal Effect of Constituent Income", 
     main="")

### Add a horizontal line for ``zero"
abline(h=0, lty=3)

### Add your lines using the ``curve" function
curve(margEff, from = xlb, to=xup, lwd=2, add=T)

##### Calculate confidence intervals
#### Get the variance-covariance matrix
covmat <- vcov(mle2)
covmat

#### Get variances/covariances
vcovBeta1 <- covmat["incomecenter","incomecenter"]
vcovBeta3 <- covmat["I(incomecenter * (hhincenter/10000))","I(incomecenter * (hhincenter/10000))"]
covBeta1Beta3 <- covmat["incomecenter","I(incomecenter * (hhincenter/10000))"]

### Make a function that calculates variance of marginal effect given Z
margEffVar <- function(z){
  return(vcovBeta1 + (z^2)*vcovBeta3 + 2*z*covBeta1Beta3)
}
### Make a function that calculates the lower bound of the 95% confidence interval
margEff_95CI_LB <- function(z){
  return(margEff(z) - qnorm(.975)*sqrt(margEffVar(z)))
}
### Make a function that calculates the upper bound of the 95% confidence interval
margEff_95CI_UB <- function(z){
  return(margEff(z) + qnorm(.975)*sqrt(margEffVar(z)))
}

##### Add the upper and lower CI curves to your plot
# Lower
curve(margEff_95CI_LB, from=xlb, to=xup, lty=2, add=T)
# Upper
curve(margEff_95CI_UB, from=xlb, to=xup, lty=2, add=T)

#### Add a vertical line for the mean of the moderating variable
abline(v=mean(mine), lwd=2, lty=2, col="red")

#### Add a rug plot
rug(mine)

#### Make a transparent histogram
#### This is a function that turns any color semi-transparent
#### ``alpha" controls the degree of transparency
makeTransparent<-function(someColor, alpha=100){
  newColor<-col2rgb(someColor)
  apply(newColor, 2, function(curcoldata){rgb(red=curcoldata[1], green=curcoldata[2],
                                              blue=curcoldata[3],alpha=alpha, maxColorValue=255)})
}

#### Make a transparent grey color
hist_col = makeTransparent("grey")

# Tell R to add any new plots to the current one
par(new=T)
# Add a histogram with no axes and with the transparent color ``hist_col"
hist(mine, axes=F, xlab="", ylab="",main="", border=hist_col, col=hist_col)
# Tell R to stop adding new plots to the current one
par(new=F)


##################################################
#       Key Vote:District Competitiveness        #
##################################################

mle2 <- lmer(keyvote ~ incomecenter + black + latino + votereported + actcenter + knowcenter + inparty + I(hhincenter/10000) + ginicenter + rephouse + unioncenter + abspvicenter + I(incomecenter*(hhincenter/10000)) + ginicenter*incomecenter +  abspvicenter*incomecenter + unioncenter*incomecenter  + rephouse*incomecenter + (1| unique.district), data = mydata2)
summary(mle2)

beta1 <- -0.344473
beta3 <-  0.040893 

margEff <- function(z){
  return(beta1 + z*beta3)
}

mine<- na.omit(mydata2$abspvicenter)



### Figure out the lower and upper x-axis bounds for plotting
xlb <- min(mine)
xup <- max(mine)

### Create a blank plot to add your lines to
plot(x=c(), y=c(), ylim=c(-1, 1), xlim=c(xlb, xup), 
     xlab="District Competitiveness", ylab="Marginal Effect of Constituent Income", 
     main="")

### Add a horizontal line for ``zero"
abline(h=0, lty=3)

### Add your lines using the ``curve" function
curve(margEff, from = xlb, to=xup, lwd=2, add=T)

##### Calculate confidence intervals
#### Get the variance-covariance matrix
covmat <- vcov(mle2)
covmat

#### Get variances/covariances
vcovBeta1 <- covmat["incomecenter","incomecenter"]
vcovBeta3 <- covmat["incomecenter:abspvicenter","incomecenter:abspvicenter"]
covBeta1Beta3 <- covmat["incomecenter","incomecenter:abspvicenter"]

### Make a function that calculates variance of marginal effect given Z
margEffVar <- function(z){
  return(vcovBeta1 + (z^2)*vcovBeta3 + 2*z*covBeta1Beta3)
}
### Make a function that calculates the lower bound of the 95% confidence interval
margEff_95CI_LB <- function(z){
  return(margEff(z) - qnorm(.975)*sqrt(margEffVar(z)))
}
### Make a function that calculates the upper bound of the 95% confidence interval
margEff_95CI_UB <- function(z){
  return(margEff(z) + qnorm(.975)*sqrt(margEffVar(z)))
}

##### Add the upper and lower CI curves to your plot
# Lower
curve(margEff_95CI_LB, from=xlb, to=xup, lty=2, add=T)
# Upper
curve(margEff_95CI_UB, from=xlb, to=xup, lty=2, add=T)

#### Add a vertical line for the mean of the moderating variable
abline(v=mean(mine), lwd=2, lty=2, col="red")

#### Add a rug plot
rug(mine)

#### Make a transparent histogram
#### This is a function that turns any color semi-transparent
#### ``alpha" controls the degree of transparency
makeTransparent<-function(someColor, alpha=100){
  newColor<-col2rgb(someColor)
  apply(newColor, 2, function(curcoldata){rgb(red=curcoldata[1], green=curcoldata[2],
                                              blue=curcoldata[3],alpha=alpha, maxColorValue=255)})
}

#### Make a transparent grey color
hist_col = makeTransparent("grey")

par(new=T)

# Add a histogram with no axes and with the transparent color ``hist_col"
hist(mine, axes=F, xlab="", ylab="",main="", border=hist_col, col=hist_col)
# Tell R to stop adding new plots to the current one
par(new=F)



###################################################
#       Key Vote: District Unionization           #
###################################################

mle2 <- lmer(keyvote ~ incomecenter + black + latino + votereported + actcenter + knowcenter + inparty + I(hhincenter/10000) + ginicenter + rephouse + unioncenter + abspvicenter + I(incomecenter*(hhincenter/10000)) + ginicenter*incomecenter +  abspvicenter*incomecenter + unioncenter*incomecenter  + rephouse*incomecenter + (1| unique.district), data = mydata2)
summary(mle2)

beta1 <- -0.344473
beta3 <-  -2.317286 # Interaction saved using : not *

margEff <- function(z){
  return(beta1 + z*beta3)
}

mine<- na.omit(mydata2$unioncenter)



### Figure out the lower and upper x-axis bounds for plotting
xlb <- min(mine)
xup <- max(mine)

### Create a blank plot to add your lines to
plot(x=c(), y=c(), ylim=c(-1, 1), xlim=c(xlb, xup), 
     xlab="District Unionization", ylab="Marginal Effect of Constituent Income", 
     main="")

### Add a horizontal line for ``zero"
abline(h=0, lty=3)

### Add your lines using the ``curve" function
curve(margEff, from = xlb, to=xup, lwd=2, add=T)

##### Calculate confidence intervals
#### Get the variance-covariance matrix
covmat <- vcov(mle2)
covmat

#### Get variances/covariances
vcovBeta1 <- covmat["incomecenter","incomecenter"]
vcovBeta3 <- covmat["incomecenter:unioncenter","incomecenter:unioncenter"]
covBeta1Beta3 <- covmat["incomecenter","incomecenter:unioncenter"]

### Make a function that calculates variance of marginal effect given Z
margEffVar <- function(z){
  return(vcovBeta1 + (z^2)*vcovBeta3 + 2*z*covBeta1Beta3)
}
### Make a function that calculates the lower bound of the 95% confidence interval
margEff_95CI_LB <- function(z){
  return(margEff(z) - qnorm(.975)*sqrt(margEffVar(z)))
}
### Make a function that calculates the upper bound of the 95% confidence interval
margEff_95CI_UB <- function(z){
  return(margEff(z) + qnorm(.975)*sqrt(margEffVar(z)))
}

##### Add the upper and lower CI curves to your plot
# Lower
curve(margEff_95CI_LB, from=xlb, to=xup, lty=2, add=T)
# Upper
curve(margEff_95CI_UB, from=xlb, to=xup, lty=2, add=T)

#### Add a vertical line for the mean of the moderating variable
abline(v=mean(mine), lwd=2, lty=2, col="red")

#### Add a rug plot
rug(mine)

#### Make a transparent histogram
#### This is a function that turns any color semi-transparent
#### ``alpha" controls the degree of transparency
makeTransparent<-function(someColor, alpha=100){
  newColor<-col2rgb(someColor)
  apply(newColor, 2, function(curcoldata){rgb(red=curcoldata[1], green=curcoldata[2],
                                              blue=curcoldata[3],alpha=alpha, maxColorValue=255)})
}

#### Make a transparent grey color
hist_col = makeTransparent("grey")

par(new=T)

# Add a histogram with no axes and with the transparent color ``hist_col"
hist(mine, axes=F, xlab="", ylab="",main="", border=hist_col, col=hist_col)
# Tell R to stop adding new plots to the current one
par(new=F)


################################
#           ADA Scores         #
################################

##########################################
#   ADA Scores: median household income  #
##########################################

par(mfrow=c(1,2))


beta1 <- -0.044243
beta3 <- -0.061635 

margEff <- function(z){
  return(beta1 + z*beta3)
}

beta3
beta1

mine<- na.omit(mydata2$hhincenter/10000)


### Figure out the lower and upper x-axis bounds for plotting
xlb <- min(mine)
xup <- max(mine)

### Create a blank plot to add your lines to
plot(x=c(), y=c(), ylim=c(-1, 1), xlim=c(xlb, xup), 
     xlab="Median Household Income", ylab="Marginal Effect of Constituent Income", 
     main="")

### Add a horizontal line for ``zero"
abline(h=0, lty=3)

### Add your lines using the ``curve" function
curve(margEff, from = xlb, to=xup, lwd=2, add=T)

##### Calculate confidence intervals
#### Get the variance-covariance matrix
covmat <- vcov(mle2)
covmat

#### Get variances/covariances
vcovBeta1 <- covmat["incomecenter","incomecenter"]
vcovBeta3 <- covmat["I(incomecenter * (hhincenter/10000))","I(incomecenter * (hhincenter/10000))"]
covBeta1Beta3 <- covmat["incomecenter","I(incomecenter * (hhincenter/10000))"]

### Make a function that calculates variance of marginal effect given Z
margEffVar <- function(z){
  return(vcovBeta1 + (z^2)*vcovBeta3 + 2*z*covBeta1Beta3)
}
### Make a function that calculates the lower bound of the 95% confidence interval
margEff_95CI_LB <- function(z){
  return(margEff(z) - qnorm(.975)*sqrt(margEffVar(z)))
}
### Make a function that calculates the upper bound of the 95% confidence interval
margEff_95CI_UB <- function(z){
  return(margEff(z) + qnorm(.975)*sqrt(margEffVar(z)))
}

##### Add the upper and lower CI curves to your plot
# Lower
curve(margEff_95CI_LB, from=xlb, to=xup, lty=2, add=T)
# Upper
curve(margEff_95CI_UB, from=xlb, to=xup, lty=2, add=T)

#### Add a vertical line for the mean of the moderating variable
abline(v=mean(mine), lwd=2, lty=2, col="red")

#### Add a rug plot
rug(mine)

#### Make a transparent histogram
#### This is a function that turns any color semi-transparent
#### ``alpha" controls the degree of transparency
makeTransparent<-function(someColor, alpha=100){
  newColor<-col2rgb(someColor)
  apply(newColor, 2, function(curcoldata){rgb(red=curcoldata[1], green=curcoldata[2],
                                              blue=curcoldata[3],alpha=alpha, maxColorValue=255)})
}

#### Make a transparent grey color
hist_col = makeTransparent("grey")

# Tell R to add any new plots to the current one
par(new=T)
# Add a histogram with no axes and with the transparent color ``hist_col"
hist(mine, axes=F, xlab="", ylab="",main="", border=hist_col, col=hist_col)
# Tell R to stop adding new plots to the current one
par(new=F)


#############################################
#   ADA Scores: District Competitiveness    #
#############################################

mle2 <- lmer(distancebasic2 ~ incomecenter + black + latino + votereported + actcenter + knowcenter + inparty + I(hhincenter/10000) + ginicenter + rephouse + unioncenter + abspvicenter + I(incomecenter*(hhincenter/10000)) + ginicenter*incomecenter +  abspvicenter*incomecenter + unioncenter*incomecenter  + rephouse*incomecenter + (1| unique.district), data = mydata2)
summary(mle2)

beta1 <- -0.044243 
beta3 <- -0.021353 

margEff <- function(z){
  return(beta1 + z*beta3)
}

mine<- na.omit(mydata2$abspvicenter)


### Figure out the lower and upper x-axis bounds for plotting
xlb <- min(mine)
xup <- max(mine)

### Create a blank plot to add your lines to
plot(x=c(), y=c(), ylim=c(-1, 1), xlim=c(xlb, xup), 
     xlab="District Competitiveness", ylab="Marginal Effect of Constituent Income", 
     main="")

### Add a horizontal line for ``zero"
abline(h=0, lty=3)

### Add your lines using the ``curve" function
curve(margEff, from = xlb, to=xup, lwd=2, add=T)

##### Calculate confidence intervals
#### Get the variance-covariance matrix
covmat <- vcov(mle2)
covmat

#### Get variances/covariances
vcovBeta1 <- covmat["incomecenter","incomecenter"]
vcovBeta3 <- covmat["incomecenter:abspvicenter","incomecenter:abspvicenter"]
covBeta1Beta3 <- covmat["incomecenter","incomecenter:abspvicenter"]

### Make a function that calculates variance of marginal effect given Z
margEffVar <- function(z){
  return(vcovBeta1 + (z^2)*vcovBeta3 + 2*z*covBeta1Beta3)
}
### Make a function that calculates the lower bound of the 95% confidence interval
margEff_95CI_LB <- function(z){
  return(margEff(z) - qnorm(.975)*sqrt(margEffVar(z)))
}
### Make a function that calculates the upper bound of the 95% confidence interval
margEff_95CI_UB <- function(z){
  return(margEff(z) + qnorm(.975)*sqrt(margEffVar(z)))
}

##### Add the upper and lower CI curves to your plot
# Lower
curve(margEff_95CI_LB, from=xlb, to=xup, lty=2, add=T)
# Upper
curve(margEff_95CI_UB, from=xlb, to=xup, lty=2, add=T)

#### Add a vertical line for the mean of the moderating variable
abline(v=mean(mine), lwd=2, lty=2, col="red")

#### Add a rug plot
rug(mine)

#### Make a transparent histogram
#### This is a function that turns any color semi-transparent
#### ``alpha" controls the degree of transparency
makeTransparent<-function(someColor, alpha=100){
  newColor<-col2rgb(someColor)
  apply(newColor, 2, function(curcoldata){rgb(red=curcoldata[1], green=curcoldata[2],
                                              blue=curcoldata[3],alpha=alpha, maxColorValue=255)})
}

#### Make a transparent grey color
hist_col = makeTransparent("grey")

par(new=T)

# Add a histogram with no axes and with the transparent color ``hist_col"
hist(mine, axes=F, xlab="", ylab="",main="", border=hist_col, col=hist_col)
# Tell R to stop adding new plots to the current one
par(new=F)


##########################
#     Simulations        #
##########################

# Do this for key vote, model 4 (mle4)

#######################################
#    Highly- Competitive Districts    #
#######################################


# rename variables
mydata2$corr.income.int <- mydata$incomecenter*(mydata$hhincenter/10000)
mydata2$district.inc <- mydata$hhincenter/10000

X <- mydata2[,c("incomecenter","black","latino","votereported","actcenter","knowcenter",
                "inparty", "district.inc", "ginicenter","rephouse","unioncenter",
                "abspvicenter","corr.income.int","IinXgi","IinXab","IinXun","IinXre")]

coeffs <- summary(mle4)$coefficients
betas <- as.matrix(coeffs[,1]) # pull out betas

varcov <- as.matrix(vcov(mle2)) # var-cov matrix for the coefficients 

set.seed(02138)
sim.betas1 <- rmvnorm(n = 10000,
                      mean = betas,
                      sigma = varcov)
dim(sim.betas1) 

Xc <- apply(cbind(1,X),2, mean, na.rm=T)

sdev <- sd(na.omit(mydata$keyvote))

# estimates at the mean 
ev.ests <- c() # Create an empty vector
for(i in 1:10000){
  miu.tilde <-Xc%*%sim.betas1[i,]
  y.ests <- rnorm(10000, miu.tilde, sdev)
  ev.ests[i] <- mean(y.ests)
}

mean(ev.ests)
quantile(ev.ests, c(.025,.975))

min(na.omit(X$abspvicenter))
min(X$district.inc)

# estimates at the mean in competitive, democratic district
ev.ests2 <- c() # Create an empty vector
for(i in 1:10000){
  Xc.new <- Xc
  Xc.new[13] <- -10.88873
  Xc.new[11] <- 0
  Xc.new[18] <- 0
  miu.tilde <- (Xc.new%*%sim.betas1[i,])
  y.ests2 <- rnorm(10000, miu.tilde, sdev)
  ev.ests2[i] <- mean(y.ests2)
}

mean(ev.ests2)

### now with the income changed
#find one standard deviation of income
sd(na.omit(mydata$incomecenter))
#3.4268

#now let's give everyone a raise
mydata2$incomecenter2 <- (mydata2$incomecenter + sd(na.omit(mydata2$incomecenter)))
mydata2$corr.int.inc2 <- mydata2$incomecenter2*(mydata2$hhincenter/10000)
mydata2$int.gini2 <- mydata2$ginicenter*mydata2$incomecenter2 
mydata2$int.abs2 <- mydata2$abspvicenter*mydata2$incomecenter2
mydata2$int.union2 <-mydata2$unioncenter*mydata2$incomecenter2
mydata2$int.rep2 <- mydata2$rephouse*mydata2$incomecenter2

X2 <- mydata2[,c("incomecenter2","black","latino","votereported","actcenter","knowcenter",
                 "inparty", "district.inc", "ginicenter","rephouse","unioncenter",
                 "abspvicenter","corr.int.inc2","int.gini2","int.abs2","int.union2","int.rep2")]

X2c <- apply(cbind(1,X2),2, mean, na.rm=T)


ev.ests3 <- c() # Create an empty vector
for(i in 1:10000){
  Xc.new <- X2c
  Xc.new[13] <- -10.88873
  Xc.new[11] <- 0
  Xc.new[18] <- 0
  miu.tilde <- (Xc.new%*%sim.betas1[i,])
  y.ests3 <- rnorm(10000, miu.tilde, sdev)
  ev.ests3[i] <- mean(y.ests3)
}

mean(ev.ests3)
quantile(ev.ests3, c(.025,.975))

## graph these

one <- density(ev.ests)
two <- density(ev.ests2)
three <- density(ev.ests3)

plot(one, xlim=c(55,70), ylim=c(0,1.2), main="",
     xlab= "Key Vote Representation (%)")
polygon(one, col=rgb(0.8,0.8,0.8,0.5))
polygon(two,col=rgb(1,0,0,0.5))
polygon(three,col=rgb(0,0,1,0.5))
#polygon(three, col=rgb(0,0,1,0,5))
abline(v=mean(ev.ests),
       col="darkgrey",
       lwd=1, lty=3)
abline(v = mean(ev.ests2),
       col="darkgrey",
       lwd=1, lty=3)
abline(v = mean(ev.ests3),
       col="darkgrey",
       lwd=1, lty=3)
temp <- locator(1)
text(temp,"Average District", cex=0.8, font=2)
temp1 <- locator(1)
text(temp1,"Avg.\nConstituent", cex=0.7)
temp2 <- locator(1)
text(temp2,"Avg.\nConsituent", cex=0.7)
temp3 <- locator(1)
text(temp3,"Wealthy\nConstiuent", cex=0.7)
temp4 <- locator(1)
text(temp4, "Democratic Competitive\nDistricts", cex=0.8, font=2)


#####################################
#       Low-Income Districts        #
#####################################

#### find top 10% and bottom 10% of districts (based off median income)
quantile(na.omit(mydata$distmedinc), c(.10,.90))
# upper = 88589
# lower = 48842
lowinc <- subset(mydata, distmedinc<=48842)
summary(lowinc$district.inc)
max(na.omit(lowinc$district.inc))
# -1.057553 <- this is what we're going to use as the value of interest

ev.ests4 <- c()
for(i in 1:10000){
  Xc.new <- Xc
  Xc.new[9] <- -1.057553
  Xc.new[11] <- 0
  Xc.new[18] <- 0
  miu.tilde <- (Xc.new%*%sim.betas1[i,])
  y.ests4 <- rnorm(10000, miu.tilde, sdev)
  ev.ests4[i] <- mean(y.ests4)
}

mean(ev.ests4)

#now with income bump

ev.ests5 <- c()
for(i in 1:10000){
  Xc.new <- X2c
  Xc.new[9] <- -1.057553
  Xc.new[11] <- 0
  Xc.new[18] <- 0
  miu.tilde <- (Xc.new%*%sim.betas1[i,])
  y.ests5 <- rnorm(10000, miu.tilde, sdev)
  ev.ests5[i] <- mean(y.ests5)
}

mean(ev.ests5)
quantile(ev.ests5, c(.025,.975))

## graph these

four <- density(ev.ests4)
five <- density(ev.ests5)

plot(one, xlim=c(55,70), ylim=c(0,1.2), main="",
     xlab= "Key Vote Representation (%)")
polygon(one, col=rgb(0.8,0.8,0.8,0.5))
polygon(four,col=rgb(1,0,0,0.5))
polygon(five,col=rgb(0,0,1,0.5))
#polygon(three, col=rgb(0,0,1,0,5))
abline(v=mean(ev.ests),
       col="darkgrey",
       lwd=1, lty=3)
abline(v = mean(ev.ests4),
       col="darkgrey",
       lwd=1, lty=3)
abline(v = mean(ev.ests5),
       col="darkgrey",
       lwd=1, lty=3)
temp <- locator(1)
text(temp,"Average District", cex=0.8, font=2)
temp1 <- locator(1)
text(temp1,"Avg. Constituent", cex=0.7)
temp2 <- locator(1)
text(temp2,"Avg. Consituent", cex=0.7)
temp3 <- locator(1)
text(temp3,"Wealthy Constiuent", cex=0.7)
temp4 <- locator(1)
text(temp4, "Democratic Low-Income\nDistricts", cex=0.8, font=2)


#####################################
#       High-Union Districts        #
#####################################

#### find top 10% and bottom 10% of districts (based off union presence)
quantile(na.omit(mydata$union), c(.10,.90))
# upper = 0.390625 
# lower = 0.140000
highunion <- subset(mydata, union>=0.390625)
summary(highunion$unioncenter)
max(na.omit(highunion$unioncenter))
# 0.3384947 <- this is what we're going to use as the value of interest

ev.ests6 <- c()
for(i in 1:10000){
  Xc.new <- Xc
  Xc.new[12] <- 0.3384947
  Xc.new[11] <- 0
  Xc.new[18] <- 0
  miu.tilde <- (Xc.new%*%sim.betas1[i,])
  y.ests6 <- rnorm(10000, miu.tilde, sdev)
  ev.ests6[i] <- mean(y.ests6)
}

mean(ev.ests6)

#now with income bump

ev.ests7 <- c()
for(i in 1:10000){
  Xc.new <- X2c
  Xc.new[12] <- 0.3384947
  Xc.new[11] <- 0
  Xc.new[18] <- 0
  miu.tilde <- (Xc.new%*%sim.betas1[i,])
  y.ests7 <- rnorm(10000, miu.tilde, sdev)
  ev.ests7[i] <- mean(y.ests7)
}

mean(ev.ests7)
quantile(ev.ests7, c(.025,.975))

## graph these

six <- density(ev.ests6)
seven <- density(ev.ests7)

plot(one, xlim=c(55,70), ylim=c(0,1.2), main="",
     xlab= "Key Vote Representation (%)")
polygon(one, col=rgb(0.8,0.8,0.8,0.5))
polygon(six,col=rgb(1,0,0,0.5))
polygon(seven,col=rgb(0,0,1,0.5))
#polygon(three, col=rgb(0,0,1,0,5))
abline(v=mean(ev.ests),
       col="darkgrey",
       lwd=1, lty=3)
abline(v = mean(ev.ests6),
       col="darkgrey",
       lwd=1, lty=3)
abline(v = mean(ev.ests7),
       col="darkgrey",
       lwd=1, lty=3)
temp <- locator(1)
text(temp,"Average District", cex=0.8, font=2)
temp1 <- locator(1)
text(temp1,"Avg.\nConstituent", cex=0.7)
temp2 <- locator(1)
text(temp2,"Avg. Consituent", cex=0.7)
temp3 <- locator(1)
text(temp3,"Wealthy Constiuent", cex=0.7)
temp4 <- locator(1)
text(temp4, "Democratic High-Union\nDistricts", cex=0.8, font=2)







